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H We introduce a new method to simulate the physics of rare events. The method, an extension of the Temperature Accelerated 

o 

Molecular Dynamics, comes in use when the collective variables introduced to characterize the rare events are either non- 
analytical or so complex that computing their derivative is not practical. We illustrate the functioning of the method by studying 
[T | the homogeneous crystallization in a sample of Lennard- Jones particles. The process is studied by introducing a new collective 
t-H variable that we call Effective Nucleus Size fA£ We have computed the free energy barriers and the size of critical nucleus, 
i- 1 which result in agreement with data available in literature. We have also performed simulations in the liquid domain of the phase 
diagram. We found a free energy curve monotonically growing with the nucleus size, consistent with the liquid domain. 

Oh 
S 

O 1 Introduction 

. rH Over the last two decades several new methods have been introduced to sample the free energy surface as a function of a set of 

£>^> collective variablesEQH These methods have been applied to many challenging problems in chemistry,^ biology^ and material 
^^^^ 

O j science.^ All these methods consist of an extended set of equations of motion coupling the dynamics of the atoms with that of 
a set of appropriate collective variables. In particular, the dynamics of the atoms is biased by a term that forces them to be in 
J> configurations compatible with the current realization of the collective variables. More explicitly, the coupling term is a function 
of the difference between the current value of the collective variable 0(x) and that of the corresponding additional dynamical 
variable z- Often the coupling is of the quadratic form (k/2) (Q(x) — z) 2 , where k is the coupling parameter. The difference among 
these methods is in how the collective variables are forced to move out of metastable states. In the Temperature Accelerated 
Molecular Dynamics^ (TAMD), of which the present method can be considered a direct extension, the variables associated 

T ^ with the collective variables are evolved at an artificially high temperature T. Since the time required to overcome free energy 
barriers is roughly proportional to exp[AF /ksT], where AF is the magnitude of barrier and ks is the Boltzmann constant, to 

^ higher T corresponds a shorter characteristic time to overcome them. By an adequate choice of T it is possible to make this 
time compatible with the maximum time achievable in atomistic simulations. A similar approach, but without the introduction 
of collective variables, has been investigated by VandeVondele and Rothlisberger 9 and by Rosso et al.E^In Metadynamics^El 
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the z are forced to visit new states by biasing the free energy in the regions already visited by these auxiliary variables. 

In all these methods the atoms are evolved by molecular dynamics and this requires the calculation of the biasing force 
—k(d(x) — z)V x 0(x). This fact limits their application to problems described in terms of collective variables which are analytical 
with respect to atomic positions. In fact, in high dimension (a very typical case), the numerical calculation of the gradient V x 0(x) 
would be computationally really challenging. However, interesting cases of collective variables which are either non-analytical 
or for which the calculation V x 0(x) is either too complex or computationally too expensive do exist. An example of the first 
case is the "ring-size" collective variable, which is used in the study of formation of clathrates (gas hydrates) to distinguish this 
kind of crystals from ice[[] 13 (which can also be formed in the same conditions). Examples of the second case are, in ab initio 
simulations, quantum mechanical observables which are not function of the Hamiltonian[§] In fact, in this case some perturbation 
theory method should be applied to calculate the biasing force. ^ 

The aim of this paper is to introduce an extension of TAMD, let us call it Temperature Accelerated Monte Carlo (TAMC), 
which allows to treat these more general cases. Since TAMC is inspired by the TAMD, and is based on the same assumptions, 
we first revise TAMD (Sec. [2]), then introduce TAMC (Sec. [5]), and finally illustrate the method by showing how TAMC allows 
to study homogeneous nucleation as described by a new collective variable, introduced in this paper, which we call Effective 
Nucleus Size (ENS) (Sec.gJ. 

2 TAMD: Temperature Accelerated Molecular Dynamics method 

In TAMD we introduce the following set of coupled equations (for simplicity, we denote x and z as scalar variables but are, 
indeed, vectors of suitable dimension): 

mx — —V x Uk(x,z)-\-thermo($) 

jli'z = -V z U k (x,z)+thermo($) (1) 

where m is the physical mass, P = (ksT)~ l and thermo($) indicates that the atoms are coupled to a thermostat at p. ju is the 
inertia of z, a parameter that can be tuned so as to achieve the adiabatic separation of the x dynamics with respect to that of z 
(see below), P = (ksT ) _1 and thermo($) indicates that the auxiliary variables z are coupled to a thermostat at p. The potential 
U k (x,z) = V(x) +k/2(Q(x) -z) 2 is the sum of the physical and a biasing potential. In the limit in which z is much slower than x, 
the force acting on z (— V z £4(x,z) = k(Q(x(t)) —z)= fk(z)) can be substituted by the time-averaged force: 

f k (z) = lim - f dtk(<d(x(t))-z) = Zj-\z) f dxk(Q(x)-z)exp[-$U k (x,z)] (2) 

f In clathrates all the water molecules are part of four, five and six-member rings. In ice, no water molecule forms rings of any size. 

§Be A(x) an observable defined as the expectation value of the operator J2(r;x) over the electronic ground state \|/(r;x) of the Hamiltonian H{r\x) (A(x) =< 
\|/|.#(r;x)|\|/ >), where r and x are the electronic and atomic coordinates, respectively. If the operator R(r;x) is a function of r and x via the !H{r;x) (J^L(r;x) = 
JZ(y{(r,x))), then, following the Hellman and Feynman theorem, =< \|/| d ^x^ |v >• However, if JZ(r;x) does not depend on x via a function of 

J-((r\x), depends also on the derivative of the ground state wavefunction with respect to x: =< \|/| |\|/ > + < \|/| R(r;x) \ ^ > + < ^ \A{r\x) \\\f >. 
This implies that if J^L(r;x) is not a function of !H(r;x), the use of the observable A (x) as a collective variable in the biased dynamics will require the calculation 
of the perturbed wavefunction . 
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where Z^{z) = f dxexp[— p£4(je,z)]. The inertia ^ can be tuned so as to obtain an adequate separation of the characteristic 
evolution times of z and x. The second equality in eq. [2] stems from the assumption that, apart for the z, the remaining degrees 
of freedom of the system are ergodic. The effective force can be interpreted as the derivative of the effective potential Fk(z) = 
— P _1 ln[Zk(z) I Z], where Z = f dxexp[— pV(x)] is the canonical partition function of the real system. Since Z is z-independent 
its introduction does not affect our argument but it is necessary for the interpretation of the effective potential Fk(z) as a free 
energy. Noting that lim^ooexp[-pf (e(x) -z) 2 ]/y/2n/$k = 8(0(jc) -z), in the limit of k -> oo F k {z) -P _1 /n[Pe(z)] = F(z), 
where Pq(z) is the probability density function that the system is in the state 0(jc) = z. In other words, for k sufficiently large, z is 
a set of random variables moving on the free energy surface and it is therefore distributed according to 

P 9 (z)=exp[-PF(z)] (3) 

Since the time average in eq. [5] is taken over the dynamics os the x's thermalized at p, the z is distributed according to the free 
energy F(z) at the physical temperature. However, since in eq. [3] the free energy F(z) is multiplied by P << p, the sampling of 
the unlikely regions is enhanced and the z can quickly overcomes the barriers on the physical free energy surface F(z). 

3 TAMC: Temperature Accelerated Monte Carlo method 

To introduce TAMC we start by observing that the key point in TAMD is that the (slow) variable z, being driven by k(Q(x) — z), 
evolves indeed according to the effective force defined in the r.h.s. of eq. [5J which is the average of k(Q(x) — z) over the canonical 
ensemble for the biased potential Uk(x,z). In TAMD this ensemble average is computed by (adiabatically separated) molecular 
dynamics (MD). However, this ingredient is not crucial and the MD on the x's could be replaced by any adequate sampling 
technique. For example, we could replace the x's MD by Monte Carlo (MC). In TAMC we take advantage of this freedom to 
replace MD by MC, which does not require the calculation of V x 0(x) for the evolution of the x and can therefore be used in 
combination with non-analytical collective variables. Let us imagine to evolve the dynamics of the collective variables according 
to the Langevin dynamics, then the TAMC can be expressed as follows 

fjz = k(Q(x MC ) -z)—yi+ (0 (4) 

where y is the friction coefficient and r\(t) is a Gaussian process with mean and covariance (r\(t)r[(t')) = 8(f — t'). In eq. 
[4] the subscript MC on the x indicates that the force &(0(xmc) — z) is computed according to the atomic configuration evolved 
by MC, in parallel to z, according to the potential Uk(x,z). If z is slow with respect to x, in a sense that will be made precise 
below, z will evolve according to the effective force < k(Q(x) —z) >u k - Following the same argument of Maragliano and Vanden- 
Eijnden, 8 in the limit k — » 00 z is distributed according to p(z) = exp[— jlF(z)]. It remains to make unambiguous the concept 
of adiabatic separation in TAMC. In TAMD, where both x and z follow a proper dynamics, adiabatic separation means that the 
characteristic time of z (x z ), i.e. the time required for a significant displacement of z, is much longer than the characteristic time 
of x (t x ): x z » x x . In TAMC x is evolved by MC and therefore the definition of adiabatic separation mentioned before cannot 
be applied. We need to introduce a different definition. Let h be the timestep used for the numerical integration of eq. [4] then 
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n z = x z /his the number of timesteps required for a significant displacement of z. let n x be the number of MC steps required for a 
accurate sampling of the probability p(jc|z) = exp[— (3C4(x,z)]/^(z) at a given z. Then, if n z » n x , z evolves according to the 
effective force < k(Q(x) — z) >u k - We stress again that, at variance with TAMD (or, for that matter, Metadynamics or Adiabatic 
Dynamics), the evolution of the atomistic configuration in TAMC does not require the gradient of the collective variable V x Q(x). 
This method, therefore, applies also to non-analytical collective variables. 

The algorithm implementing the method described above is quite simple and, in practice, consist in a dynamics over the 
collective variables z with the force k(d(xMc) — z) computed according to the atomic configuration evolved by a Metropolis 
Monte Carlo controlled by the potential C4(x,z). This algorithm is presented schematically in Fig. [T] Starting from the atomic 
configuration x l and the value z l of the collective variable and its velocity z l , the values z* +1 / 2 at the "half time step" and z l at 
the next time step are computed according to the force k(d(x l ) — z l ). After this, the atomic position is evolved according to the 
Metropolis MC scheme: i) an atom is chosen at random, ii) a random displacement is applied on this atom so obtaining the 
configuration x*, iii) the variation in the (biased) potential energy AU^ = £4(x*,z i+1 ) — Uk{x\z l+l ) is computed, iv) the move 
is accepter (rejected) according to the usual Metropolis criterion: either with probability 1 if the AUk < or with a probability 
equal to exp[— PA£4]. The only difference with respect to the standard MC move is that in TAMC the variation of energy is the 
sum of the variation of the interatomic potential, which is computed according to the standard procedure, plus the sum of the 
biasing term k/2(Q(x) — z) 2 . This procedure is repeated for Nmc steps and the force < k(Q(x) —z) >u k is computed as the average 
of k(Q(x) — z) over the configurations generated by the MC sampling. We would like to remark that any method developed for 
standard MC can also be applied to TAMC. In particular, we can compute not only the Helmhotz free energy but also the Gibbs 
free energy by running constant pressure MC simulation, as it was done in this paper to study the crystal nucleation (see Sec. 

We now compare TAMD and TAMC from the point of view of the various parameters governing the simulation. In TAMD 
the timestep h is fixed by the timescale of the fast degrees of freedom, that is the atoms. Typically, the timestep is of the order of 
10 -15 s (1 fs). The parameter ju, the inertia of the variable z, is adjusted so as to achieve the adiabatic separation of the dynamics 
of the z and the x. Depending on the type of thermostat, other parameters must be adjusted so that atoms are adiabatically at the 
thermal equilibrium. For example, if eq.s[T]take the form of Langevin equations, the friction parameters y and Y, for the x and 
the z, need to be adjusted. If eq.s[T]take the form of Nose-Hoover equations the inertia Q and Q' of the thermostat variables need 
to be defined. These parameters control the timescale over which in eq. [2] we can assume the time average to be equivalent to the 
ensemble average. In TAMC we must chose the MD and MC parameters controlling the evolution of z and x, respectively, such 
that the adiabatic separation as described before (i.e. n z » n x ) is achieved. The MC on the x is governed by the usual parameters 
(maximum atomic displacement for standard MC) according to the usual criteria on the acceptance ratio (typically between 30 
and 60 %). The choice of these parameters determine the number of MC steps n x needed between the Langevin timesteps of the 
z. Once n x is determined, the inertia /i, thermostat parameters (e.g. y for Langevin and Q for Nose-Hoover) and the length of the 
MD timestep are set such that the adiabatic separation is achieved. 
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Fig. 1 Flowchart of the TAMC algorithm. In this chart we assume that the MD step is performed by a three-step algorithm, in particular we 
used the Vanden-Eijnden and Ciccotti algorithm for integrating Langevin equation of motion 



4 An example of application of TAMC: homogeneous nucleation in supercooled Lennard- Jones liq- 
uids. 

In this section we shall illustrate how TAMC works in practice by studying the homogeneous nucleation in a Lennard- Jones 
liquids. The nucleation is studied as a function of one collective variable monitoring the size of the crystalline nucleus (see 



Sec. 4.1 ). Since we use only one collective variable, other method could be used as well to tackle this problem. Indeed, the 
nucleation in Lennard- Jones liquids, using different collective variables, has been already studied using the Umbrella sampling 
and the Partial Path Transition Interface Sampling techniques (see Refs. [15] and [17]). 

4.1 Collective variable for nucleation 

One collective variable proposed to study homogeneous nucleation of liquids is the size $£{r\ , . . . , r#) of the largest cluster of 
solid-like particlesESEH] expressed in terms of number of particles. In the present work we found more suitable to adopt an 
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improved definition of the size of a crystalline nucleus. Our definition, and its advantages over that used in the referred works, 
are given in the following. As in previous works 1 -^H to define fA£(ri , . . . , r#) we start from the local bond-orientational order 
parameter of Steinhardt et al.E3 

Ni 
7=1 

where 7/ m (r iJ ) is the spherical harmonics of order / and m computed at the polar and azimuthal angles associated to the vector 
?ij connecting atoms i and j. The sum runs over the Nf nearest neighbors of the atom i. Here and in the following, the nearest 
neighbors of an atom, say i, are those atoms satisfying the condition |r^| < 1.5 a. For a proper choice of q\ m (i) characterizes 
the local structure around the atom i. In crystals the qi m of neighboring atoms are similar (coherent). In liquids the degree of 
coherence, in a sense made precise below, is much lower than in crystals. We will use this empirical observation to distinguish 
between crystal-like from liquid-like particles. The degree of coherence of q\ m of neighboring atoms can be measured by 

\l! m =-iqim(i)*qimU)\ 

U ~ \qi®\\qiU)\ W 
where \qi(i)\ and \qi(j)\ are defined as \qi{i) \ = ^Y} m =-i<lim(i)*qim(i). <27m(0* is the complex conjugate of q\ m (i). In fact, Qj 
can be interpreted as the normalized dot product between the vectors qi{i) and qi(j), with components qi m (j). Therefore Qj 
measures how much the two vectors are parallel. When two atoms have the same environment their q\ m are the same and Qj = 1. 
This is the case of two neighboring atoms in a crystal at K. On the contrary, when two atoms have a different environment, 
as it is typically the case in disordered systems, Qj is much smaller. At finite temperature Qj ^ 1 also in crystals. However, 
it has been empirically observed that Qj is > 0.5 in crystalline samples at finite temperature. 16 According to this, we say that 
two particles with a Qj above this threshold are "connected". We shall use throughout the paper the word "connected" in this 
strict sense. Neighboring particles with an alike environment (Qj > 0.5) are found also in liquids, but in crystals the number 
of connected particles around a given atom (typically > 7) is higher than in liquids (typically below 4 — 5). So, we can use the 
number of connected particles as a parameter to distinguish between liquid-like and crystal-like particles. 

Once we have identified crystal-like particles, we can search for clusters of such particles in the sample. In Refs. [15] and [17] 
a crystal-like particle constitute, or is assigned to, a crystal nucleus if it is within a given distance (1.5 a) from another crystal-like 
particle (possibly of an already existing nucleus), whether this last particle is connected or not to it. On this, our definition of 
crystal nuclei is different. We define a nucleus as the set of crystal-like particles which are connected (and therefore close, in 
the sense defined above) with at least another particle in the cluster (see Fig. [2]). According to the definition given in Refs. [15] 

f A first aspect that must be taken into account in the choice of / is that spherical harmonics corresponding to odd / are anti-symmetric under inversion (F/ m (?) = 
— Yim(— ?))• in many crystals the nearest neighbors of an atom i have an equilibrium configuration such that if one of them is connected to i by the vector ry, 
there is another one connected by — (all the simple lattices - simple cubic, FCC, BCC, etc - have this property). This implies that the qi m 's with odd / of 
these systems are all 0. Then, for qi m to be useful in distinguishing between different local environments, / must be even. Steinhardt et al. 19 for bulk and cluster 
systems at their equilibrium configuration and ten Wolde et alM^for nucleation have shown that the optimal / to use for qi m and derived observables (introduced 
to distinguish between systems of different symmetry) are / = 4 or / = 6. In this paper we use / = 6 to identify connected particles, as this is consistent with the 
/ used in Ref. [17]. That choice allows a direct comparison with previous results. 
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Fig. 2 Schematic representation of atoms forming clusters and their connectivity. The clusters 'A' and 'B' are two separated clusters as no one 
atom in one cluster is connected with any atom in the other, even though there are atoms of one cluster which are nearest neighbors of atoms 
of the other cluster. This is because the orientation of atoms in one cluster are different from the one in the other and upon growing they will 
not form a single crystal. 

and [17] the nucleus might include also particles which are not connected, i.e. particles which have a different environment. On 
the contrary, our definition of clusters allows to distinguish between particles which are indeed members of the same nucleus 
from particles which might be member of neighboring nuclei. The two different nuclei, if growing, will eventually form a grain 
boundary, like the nuclei 'A' and 'B' of Fig. [2] We expect that by adopting our definition of crystal nucleus the size of critical 
nuclei observed in Ref. [15] and [17] would be reduced in size by ~ 10 — 50 atoms, which is the typical size of subcritical crystal 
nuclei usually found in Lennard-Lones liquids (see top panel of Fig. [5]) and that could be close in space to the largest connected 
crystalline cluster. Indeed, the fact that the definition of Ref. [17] might include in the crystalline nucleus also particles that 
are not connected to it could account for the broad shape of the committor observed in this paper. In fact, as explained before, 
in this case clusters of significantly different sizes, both under and super-critical, could be assigned to the critical nucleus. In 
practice, we identify these clusters by using methods of the Graph theory. In our approach the crystal-like particles are the nodes 
and the connections between them are the edges of a graph. We then identify clusters in the graph by means of the Deep First 
Search method. 

m xhe 

Deep First Search method consists in searching the graph for connected nodes starting from a root node, 
and exploring as far as possible along the edges of the graph. Once all the particles directly or indirectly connected to the root 
are identified the search for the members of the present cluster is completed. The search for the particle belonging to another 
cluster is then started by defining a new root among the crystal-like particles not yet assigned to any nucleus. This search is 
repeated until all the crystal-like particles have been assigned to a cluster. The collective variable fA£(n , • • • , r#) could be defined 
as the number of particles forming the largest cluster. It appears clear from the description above that this collective variable is 
non-analytical as it does depend on atomic positions r but not through an explicit formula, rather through a search procedure. 
Indeed both definition of nucleus size, the present one and the one given in Refs. [15] and [17], are non- analytical. 

The definition of the collective variable fAt(n , . . . , r#) given above is not very efficient to use with TAMC since it is discrete 
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Fig. 3 Nucleus and buffer particles in a sub-critical nucleus. For sake of readability, nucleus particles are bigger (and green online) and buffer 
particle are smaller (and red online). The sum of eq. j7Jruns over both sets of atoms. 

and therefore the associated potential (k/2) (fA£(n , • • • , r#) — z) 2 biases the simulation only when the MC move changes its value. 
When this does not occur, the MC is unbiased and there is no acceleration on the sampling. Let us illustrate this problem by 
an example. Let us assume that fA£(n , • • • , r^) is lower than z. This means that a MC move increasing fA£(n , • • • , r#) should be 
favored by the biasing potential. An MC move increases fA£ (r\ , . . . , r#) if it connects a crystal-like particle not yet member of the 
nucleus to another particle which is in the nucleus or if it transforms a non crystal-like particle connected to the nucleus into a 
crystal-like one. However, both processes might require several steps. For example, the first might occur through a series of MC 
moves which gradually increase the degree of crystallinity of a non crystal-like particle (the meaning of degree of crystallinity, 
made precise below, for the time being must be understood as the number of particles connected to the present one). However, 
as far as this particle does not become crystal-like, all these moves are unbiased, which slows down the convergence of the 
MC procedure. Moreover, since the definition of connectivity and crystallinity of a particle is somewhat arbitrary, as it depends 
on the value of the lower bound of Qj for which we consider two particles connected, the current value of the nucleus size is 
strongly dependent on this arbitrary choice. We solved the problem of biasing all the MC moves and alleviated the problem of 
the arbitrariness of fA£(n , . . . , r#) by introducing a continuous version of the collective variable described above. This continuous 
version of the collective variable we have called Effective Nucleus Size. It is defined by the following equation 

fA£(n,...,r^)= £ w c (i)w n (i) (7) 

iticluster+buffer 

where the sum runs over the members of the largest cluster as defined above plus the atoms laying within a 1.5 a thick buffer 
around the cluster (see Fig. [3} and w c {i) and w n (i) are weights accounting for the degree of connection of the particle i to the 
cluster and its degree of crystallinity, respectively. The buffer particles are those particles which do not belong to the nucleus and 
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which are within a distance of 1.5 a from any nucleus particle. The so defined fA£(n , • • • , r^) is still non-analytical as the sum in 
eq. ([7} can be defined only after the atoms in the largest cluster are identified through the procedure explained before. 




Qmax 

1[ ' 1 ' 1 i f • f • f 

0.8 - • 



•3 0.6 

0.2 - • 

()• • 4 * 1 1 ■ 1 1 ■ 

2 4 6 8 10 12 

n(i) 

Fig. 4 w c and w w weights as a function of the degree of connection C™ ax and the degree of crystallinity n (see text). 



In order to give an explicit expression for w c {i) we need to analyze what are the expected properties of this parameter. w c {i) 
must tend to (to 1) when the particle i is loosely (strongly) connected to the cluster, and must smoothly go from one limit to the 
other in between. In other words, it should be a sigmoid with respect to the degree of connection of the particle i to the cluster. 
The parameter accounting for the degree of connection in our modeling of nucleation is C™ ax = maxy(Qy), where j is a particle 
belonging to the cluster. A reasonable function w c {i) of C™ ax with the properties mentioned above can be obtained from the 
Fermi function 

^' ) = 1 -expMcr x -AO] + l (8) 
where ju c is the parameter controlling the value of C™ ax at which the function switches from low values to high values and X c is 

a parameter controlling the smoothness of the switching. We set ju c = 0.5. X c was set to 15 so that w c (i) is very low (< 0.05) 

for C™ ax < 0.3 and close to 1 (> 0.95) for C™ ax > 0.7 (see the top panel of Fig. [4|. Similarly, w n (i) is defined according to the 

following expression 

w »(0 = 1 r . 7 /•! vTTT ( 9 ) 

explA n (n(i)-ju n )\ + l 

where n{i) is the number of neighbors connected to the particle i. n(i), as announced above, measures the degree of crys- 
tallinity of a particle. The ju n and X n are set to 5 and 1.5, respectively. With this choice, w n (i) is close to 1 (0.95) for n(i) = 1 (see 
the bottom panel of Fig. [4]), which is consistent with the lower bound used to identify crystal-like particles for the identification 
of nucleus. 

The continuous formulation of the collective variable fA£(n , • • • , tn) given above solves the problems of the original definition. 
In fact, with this new definition, essentially any MC move is biased. In order to illustrate this, let us imagine that the current 



9 



value of 0\C (r\ , . . . , r#) is lower that the target value z. Any MC move involving atoms belonging to the largest cluster or to the 
buffer region (see Fig. [3} which increases the connection of this particle with the cluster or increase its degree of crystallinity 
will increases 9\C( r i > • • • ? r v) and therefore will be favored as they reduce the biasing potential, even though the moved particle 
is not yet crystal-like according to the definition given before or it is not yet connected to the nucleus. On the contrary, any 
MC move reducing the connection of the particle to the cluster or its degree of crystallinity will be disfavored. In other words, 
since the change of the collective variable is continuous (indeed almost continuous, since it is continuous for the term depending 
on the connectivity but not for the one depending on the number of neighbors, which is function of the discrete variable n(i)) 
the possibility of accepting MC moves that bring the system closer to the target value is higher than with the original discrete 
formulation of the collective variable $£(r\ , . . . , r#). This definitely improves the efficiency of the sampling. 

4.2 Results and Discussion 




le+05 2e+05 3e+05 4e+05 5e+05 
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Fig. 5 Timeline of the collective variable fAC(ri , . . . , r#) obtained from standard MD simulations (top panel) and by TAMC (bottom panel). 

We ran unbiased and TAMC simulations of a liquid sample of 3456 Lennard- Jones particles. The pressure and the temper- 
ature were kept fix at P = 5.6 T = 0.92, i.e. we performed isobaric-isothermal MD/MC simulations (pressure and temperature 
are reported in Lennard- Jones units). These conditions correspond to a 17% degree of supercooling, i.e. at this pressure the 
temperature is 17% lower than the melting temperature (we used the Hansen and Verlet 21 data to estimate the melting tempera- 
ture). This pressure and temperature are in the range used in other computational investigations of nucleation in Lennard- Jones 
sy stems. EHED The liquid sample is obtained by melting a Body Centered Cubic crystal at high temperature (T = 50) and then 
cooling it down to T=0.92 very slowly. 

We performed unbiased MC and MD simulations of the liquid and monitored the fA£(n , . . . , r#) (see the top panel of Fig. [5] 
where only the MD results are shown). In these simulations, on a trajectory of 55000 MD timesteps (or 2 x 10 6 steps for MC) 
we only observed the $£(r\ , . . . , r#) to fluctuate around 20. Nuclei of this size are under-critical, that is their size is smaller than 
the size corresponding to the maximum of the free energy versus !J\[ curve, which Moroni et al.ES, using their definition of fA£ 
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have found to be ~ 250 for a Lennard- Jones liquid in similar conditions (25% degree of supercooling). This result confirms that 
nucleation is a rare event and that relevant information on this process, such as the free energy barrier and the critical size of the 
nuclei, cannot be obtained by brute force simulations in the conditions of "moderate" supercooling. 
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Fig. 6 Free energy vs. fA£ curve at P=5.6 and T=0.92. The nucleus critical size is also reported. 

Starting from the same liquid sample, we run a TAMC simulation using fA£(n, . . . , r#) as collective variable. In these simu- 
lations we evolved the dynamical variable fA£ associated to fA£(n , • • • , r#) (see Eq. ^ every 3456 MC steps on the nuclei. This 
roughly corresponds to make one fA£ move every one move of all the atom. As a first remark, we notice that with TAMC we are 
able to explore a wide range of the collective variable space, as shown in the bottom panel of Fig. [5] In particular, we notice that 
TAMC allows to explore the under-critical as well as the post-critical domain of the nucleus size. 

A comment is in order about the results shown in Fig. [5] we notice that the system fluctuates between what appear to be 
metastable states in the post-critical domain. It might appear surprising that post-critical relatively stable nuclei do exist of a 
size smaller than the total number of atoms. We propose two possible explanations of origin of this phenomenon. On the one 
hand, a further growth of clusters of effective nucleus size fA£ = 400 — 600, containing of the order of one half of the particles in 
the sample, requires a proper orientation of the nucleus with respect the simulation box otherwise the mismatch might prevent 
the formation of a perfect crystal containing all the particles in the sample. The re-orientation of a nucleus of this size is a slow 
process and therefore, in absence of any acceleration, do not take place over the time scale of a simulation. Another possible 
explanation is connected to the absence of additional collective variables controlling the global level of order in the growing 
cluster. The collective variable $£(r\ , . . . , r#) alone cannot control this phenomenon as it depends only on the value Qj of the 
nearest neighbor particles and therefore it controls only the local level of ordering. Moroni et al.E^ have tried to solve this 

1 /2 

problem by adding to their 9^{r\ , . . . , r#) the collective variable Q£ l =4%/ 13 (j^ = _ 6 \ Q^ m , where Q$ m is the average of 
the local bond-orientational order parameter q$ m (i) computed over the atoms belonging to the nucleus. However, Moroni et al., 
and Gasser et al., 22 who have performed experimental work on the nucleation of colloidal particles, have shown that the free 
energy barrier and the critical nucleus size can be accurately computed (measured) taking only into account the fA£(n , . . . , r#) 
collective variable. 
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We now turn to the reconstruction of the free energy curve vs. fA£ Instead of running one long TAMC simulation we took 
advantage of parallel computers by running 32 independent TAMC trajectories. The 32 simulations were started from four 
configurations extracted from the TAMC trajectory shown in Fig. [5] Two of these four configurations correspond to under- 
critical nuclei (fA£ =100 and fA£ = 200) and two to post-critical nuclei (fA£ = 300 and fA£ = 400). For each configuration we 
started eight TAMC trajectories with different initial random values of t These initial velocities were sampled from a Maxwell- 
Boltzmann distribution at T = 30 (the nucleation barrier as estimated by Moroni et al. at T = 0.83 is AF* = 25ksT). Each 
trajectory is evolved for 150000 collective variable steps. In Fig. [6] we report in logarithmic scale the histogram of P(9Q 
(— log[P(fA0] = P^X^AQ) vs - y£ computed along the TAMC trajectory. On the basis of this curve we can measure the nucleation 
free energy barrier and estimate the critical nucleus size. The free energy barrier we measure is AF* ~ 30.5ksT, slightly higher 
than the free energy barrier computed by Moroni et al. 

ESl for 

a system at T=0.83 (corresponding to a 25% degree of super-cooling) 
by the Partial Path Transition Interface Sampling method 23 (PPTIS). It is intuitive to expect that to an higher super-cooling 
corresponds a lower free energy barrier, being the liquid state less stable in this condition. The decrease of the free energy barrier 
with the degree of super-cooling is indeed also a result of the classical theory of nucleation. ISlThis hypothesis is further supported 
by the fact that crystallization in highly supercooled systems can be studied by "brute force" simulations,^ which means that the 
free energy barrier, and therefore the corresponding timescale, decreases with the degree of supercooling. 

In spite of the difference in the free energy barrier, the critical size as estimated from our TAMC simulation is in good 
agreement with the one reported by Moroni et al.^The critical size, defined as the !J\[ corresponding to the maximum of the free 
energy, is ~ 230 in our simulations while it is 243 in Ref. [15]. Both, our and Moroni et al. results are strongly different from 
the critical size reported by ten Wolde et al.E^(fA£* = 642). However, in Ref. [17] the collective variable used for studying the 
nucleation is a different one (it is the bond-orientational order parameter of Steinhard et al. 19 of the entire sample) and the 
critical nucleus size is computed as the size of the maximum among the largest nuclei of crystal-like particles at the value of 
corresponding to the maximum of the free energy. 

We also analyzed the structure of sub-critical and post-critical nuclei. We found that sub-critical (small) nuclei are, globally, 
disordered. An example of such a nucleus is shown in Fig. [TjA. This nucleus does not show any degree of global ordering. 
Rather, it looks liquid-like. This observation is partly in line with the results of Ref. [15] and Ref. [17], where it is reported 
that small nuclei are mainly liquid- and bcc-like. On the contrary, post-critical nuclei show a more ordered structure, with an 
fcc/hcp-like core. In Fig. [TjB it is shown the ordered core of a post-critical nucleus, where the close-packed hexagonal crystal 
plane is clearly visible. This result is also consistent with recent experimental results 22 on the nucleation of colloidal particles, 
which are characterized by post-critical nuclei with a rough surface and an ordered face centered cubic or hexagonal close-packed 
core. 

Finally, we also run simulations at P = 5.6 and T = 1 .2. According to the phase diagram determined by Hansen and Verlet by 
MC simulations, 21 at this pressure and temperature the most stable phase is the liquid phase. We therefore expect the free energy 
curve to be minimal at low fA£ and monotonically increasing with the nucleus size. Indeed, in our simulations in these condition 
the free energy is minimal at low fA£, even though it seems to reach a plateau at large values. This is most likely due to the 
insufficient statistics we get for event of very low probability. However, our results confirm that at this pressure and temperature 
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(A) 



(B) 



Fig. 7 (Color online) Snapshots of an under-critical nucleus (A) and the ordered core of a post-critical nucleus (B). In panel (B) the network 
of bonds is superimposed to the atomistic structure to emphasize the ordered shape of the structure. 

the most stable phase is the liquid phase. 




Fig. 8 Free energy vs. fA£ curve at P=5.6 and T=1.2. 

5 Conclusions 

In this paper we have introduced a new method for simulating rare events described by non-analytical collective variables of 
atomic positions. This kind of variables is intrinsic to ab-initio simulations, were the collective variables are the expectation 
value of the associated operators computed over the wavefunction corresponding to a given atomic configuration. However 
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important cases of non-analytical variables appear also in classical simulations, as our illustration with the nucleation case has 
shown. 

To illustrate the functioning of the method we studied the homogeneous crystallization in a sample of Lennard- Jones particles. 
The process has been studied using the new collective variable, Effective Nucleus Size fA£ , introduced in this paper. Our results 
at P = 5.6 and T = 0.92 are in agreement with previous simulations. Moreover, simulations at a pressure and temperature in 
the liquid domain (P = 5.6 and T = 1.2) found a free energy curve growing with the nucleus size, in agreement with what we 
expect in the liquid domain. Our conclusion is that the method is ready for challenging applications. Work is in progress in this 
direction. 
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